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The intrinsic flux of very high-energy (VHE, energy ^100 GeV) 7-rays from extragalactic sources 
is attenuated due to pair production in the interaction with photons of the extragalactic background 
light (EBL). Depending on the distance of the source, the Universe should be opaque to VHE 
photons above a certain energy However, indications exist that the Universe is more transparent 
than previously thought. A recent statistical analysis of a large sample of VHE spectra shows that 
the correction for absorption with current EBL models is too strong for the data points with the 
highest attenuation. An explanation might be the oscillation of VHE photons into hypothetical 
axionlike particles (ALPs) in ambient magnetic fields. This mechanism would decrease the opacity, 
as ALPs propagate unimpeded over cosmological distances. 

Here, a large sample of VHE 7-ray spectra obtained with imaging air Cherenkov telescopes is 
used to set, for the first time, lower limits on the photon- ALP coupling constant over a large 
range of ALP masses. The conversion in different magnetic field configurations, including intra- 
cluster and intergalactic magnetic fields together with the magnetic field of the Milky Way, is 
investigated taking into account the energy dependence of the oscillations. For optimistic scenarios 
of the intervening magnetic fields, a lower limit on of the order of 10~^^GeV~^ is obtained, 
whereas more conservative model assumptions result in Qa-^ ^ 2 x 10~^^GeV~^. The latter value is 
within reach of future dedicated ALP searches. 

PACS numbers: 14.80.Va, 98.54.Cm, 95.85.Pw, 98.70.Vc 



|manuel.meyer@physik.uni-hamburg.de| 



2 



I. INTRODUCTION 

Very high-energy (VHE, energy ^100 GeV) 7-rays from extragalactic sources interact with photons of the extra- 
galactic background hght (EBL) and produce e+e~ pairs [1-3 . As a result, the intrinsic 7-ray flux is exponentially 
suppressed with the optical depth, r(z, E), that increases with both the redshift z of the source and energy E. Above a 
certain energy, the Universe should thus become opaque to 7-rays from sufficiently distant sources. The optical depth 
also depends on the photon number density of the EBL whose exact level remains unknown, as direct observations 
are difficult due to contamination by foreground emission [4 . The EBL ranges from ultraviolet (UV) / optical to 
far-infrared wavelengths and originates from the starlight integrated over all epochs and the starlight absorbed and 
reemitted by dust in galaxies (see e.g. Refs. [5l [6] for a review). Nevertheless, firm lower limits in the optical and 
infrared exist which can be derived from galaxy number counts [T, , and recent EBL models show a good agreement 
in the overall shape of the spectrum (e.g. Refs. |9HTT]). In addition, VHE-7-ray spectra from extragalactic sources 
can be used to place upper limits on the EBL density if certain assumptions about the intrinsic spectra are made (see 
e.g. Refs. [T2[[T3; for recent upper limits). 

However, indications have been found that the Universe is more transparent to VHE 7-rays than the prediction 
of EBL models suggests (e.g. Refs. p!4HT7] ). The authors of Ref. [18 introduced a nonparametric test to quantify 
the accordance between EBL model predictions and VHE observations. They tested the EBL model of Ref. JO] 
(henceforth the KD model) which is designed to produce a minimal attenuation as it closely follows the galaxy 
number counts in the infrared. The authors used a statistical sample of the spectra of active galactic nuclei (AGNs) 
measured with imaging air Cherenkov telescopes (lACTs) and fitted analytical functions to the spectral points [19 . 
Subsequently, the authors tested if the residuals of data points that correspond to an optical depth r ^ 2 follow 
the expectation to scatter around a zero mean. They found, though, that the mean is larger than zero at a ~ 
4(7 significance level, indicating an overcorrection of the data with the KD model. A possible solution to this 
observed pair production anomaly might be the conversion of photons into axionlike particles ( ALPs) [16, J/7 . ALPs 
are pseudo Nambu-Goldstone bosons that are created if additional global symmetries to the standard model are 
spontaneously broken (see e.g. Ref. [20] for a review). Such fields are a common prediction in compactified string 
theories [21 . Their phenomenology is closely related to that of axions which solve the strong CP problem in QCD 
[22H24] . Most importantly in the present context, ALPs share the same coupling to photons as axions, characterized 
by the Lagrangian 



Ca^ = -^9a^ E^uF^^a = ga^E-Ba, (1) 

where E^^ is the electromagnetic field tensor (with electric and magnetic fields E and B, respectively), E^^ is its 
dual, a is the ALP field strength, and is the photon- ALP coupling strength which has the dimension (Energy)"^. 
In contrast to the axion, the ALP mass is unrelated to the coupling strength. 

The effect of photon- ALP mixing on VHE 7-ray spectra in different magnetic field settings has been extensively 
discussed in the literature. The case of a conversion in an intergalactic magnetic field (IGMF) was addressed by, e.g., 
Refs. p!6| fT7[ [25H27] . The authors of, e.g., Refs. [28l [29] included the magnetic fields in and around the source, and 
in Ref. [30] the backconversion of ALPs into photons in the galactic magnetic field (GMF) of the Milky Way was 
studied. Recently, the photon- ALP mixing for sources located inside galaxy clusters and the reconversion in the GMF 
was also investigated [31]. Usually, previous studies used fixed values for the ALP mass and coupling close to current 
experimental bounds in order to maximize the effect on the 7-ray spectra. 

In this article, the goal is to determine the preferred region in the {nia^ga-f) parameter space that minimizes the 
tension between data and model predictions found in Ref. [18 . This allows, for the first time, to place a lower limit 
on the photon- ALP coupling to explain the observed transparency of the Universe for VHE 7-rays. Four different 
scenarios for the intervening magnetic field will be considered, including mixing in the IGMF, the intracluster magnetic 
field (IGMF), and the GMF of the Milky Way. In two cases, the parameters of the IGMF and IGMF wih be chosen 
as optimistically as possible in order to derive lower limits on ga^. Additionally, a more conservative choice of 5-field 
model parameters will be investigated. Furthermore, two different EBL model realizations will be studied. 

The paper is organized as follows: In Sec. |n| the theory of photon- ALP conversions is briefly reviewed. In Sec. 
|IIH the different magnetic field configurations and the scenarios that are analyzed are summarized. Subsequently, the 
method of constraining the ALP parameter space is laid out in Sec. |IV| before presenting the actual results in Sec. |V] 



and concluding in Sec. VI 
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II. PHOTON-ALP CONVERSION IN MAGNETIC FIELDS 



The photon- ALP interaction is described by the Lagrangian 

£ = Ha-i + >CeH + f^a (2) 

where La is given in Eq. 0, >Ceh is the effective Euler-Heisenberg Lagrangian accounting for one-loop corrections in 
the photon propagator (see, e.g., Ref. [32 ) and Ca includes the kinetic and mass term of the ALP, 

Ca = ^d^ad^a - ]^mla^ . (3) 

For a monochromatic photon / ALP beam of energy E propagating along the x^, axis in a cold plasma with a 
homogeneous magnetic field, it can be shown that C leads to the following Schrodinger-like equation of motion [32]: 

where Ai{x^) and A2{x^) describe the linear photon polarization amplitudes along xi and X3, respectively, and a{x^) 
denotes the ALP field strength. Let denote the magnetic field transverse to the beam propagation direction. If 
one chooses Bt to lie only along the X2 direction, the mixing matrix M.q can be written as 



>to = A,^ , (5) 




where a mixture of the photon polarization states due to Faraday rotation can be safely neglected for the energies 
considered here. The matrix elements A|| = Api + 7/2Aqed and A^ = Api + 2Aqed account for medium effects on the 
photon propagation, where Api = —uo^\/ with the plasma frequency of the medium, uj^\. The plasma frequency is 
connected to the ambient thermal electron density nei through uj^\ = 3.69 x lO""*^"^ Y^nei/cm~^ eV. The QED vacuum 
birefringence effect is included in Aqed B'^. Furthermore, A^ = —m?^/{2E) and the photon- ALP mixing is 
induced by the off-diagonal element Aa7 = l/2ga^B±_ (see, e.g., [31 for numerical values in suitable units of the 
matrix elements). If photons are lost due to the interaction with the EBL, the elements A|| ^ are modified to include 
a complex absorption term, A||^^ ^||,^ +z/(2Ai^^P), where AJ^^p is the mean free path for photons undergoing pair 
production. Equation Q is solved with the transfer function T{x^^{)]E) with the initial condition T(0,0;^) = 1. 
Neglecting the birefringence contribution for a moment, it can be shown that the photon- ALP oscillations become 
maximal and independent of the energy E and ALP mass for an energy above the critical energy 

which is the so-called strong mixing regime. However, as the goal of this paper is to constrain the (m^, ^07) parameter 
space, it is generally not the case that the mixing occurs in this regime. 

So far, only a polarized photon beam has been considered. As of today, the polarization of VHE 7-rays cannot 
be measured and one has to consider an unpolarized photon beam and reformulate the problem in terms of density 
matrices. The general polarization matrix is given by 

P{X3) = A2{X3) (g) {Ai{x3) A2{X3) a{x3)) (7) 

V a{x3) J 

and the equation of motion takes the form of a von Neumann-hke equation, 

i^ = [P,-Mo], (8) 

which is solved by p{xs,E) = T{xs^O;E) p{0)T\xs,0; E). In the more general case in which B^ has an arbitrary 
orientation and forms an angle tjj with the X2 axis, the solution can be found via a similarity transformation 

'^cos t/j — sin t/j 0\ 

V{iIj) = I sinV^ cosV^ , (9) 
1/ 
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so that M = V{iIj)MoV'^{iIj) and the solution to the modified Eq. ^ is r{x3,0; E;iIj) = V{ilj)r{xs,0; E)V'^ {i/j). If, 
moreover, the beam path can be spht up into n domains with a constant and homogeneous magnetic field in each 
domain but a changing orientation (and strength) from one domain to the next, the complete transfer matrix is simply 
given by the product over all domains, 

n-l 

T(x3,n,^3,0;^; • • • , V^o) = Yl (^3,fc+l , ; V^/c) , (10) 

k=0 

with one mixing matrix A4k for each domain. The transition probability of observing a photon / ALP beam in the 
state pfinai after the crossing of n magnetic domains reads 

^final = Tr(pfinalT(x3,n, 2:3,0; E; ?An-l, • • • , V^o)p(^3,o)T'^ (^3,n , 2:3,0; E; l/j^-l, • • • , V^o))- (H) 

Equipped with this formula, the photon transition probability P^^ is defined as the sum of the transition prob- 
abilities from an initially unpolarized pure photon state Punpoi = l/2diag(l, 1, 0) to the final polarization states 
pii = diag(l, 0, 0) and P22 = diag(0, 1, 0): 

— Pii + P22 = Tr((pn+p22)rpunpoir^). (12) 

III. MAGNETIC FIELD CONFIGURATIONS AND SCENARIOS 

When the photon / ALP beam propagates towards Earth, it crosses different regions of plasma and magnetic field 
configurations. The following environments are considered, ordered by increasing distance from Earth: 

1. The Galactic magnetic field of the Milky Way (GMF). 

2. The intergalactic magnetic field (IGMF). 

3. The magnetic field inside a galaxy cluster (intracluster magnetic field, ICMF) in the vicinity of the emitting 
source. 

The observational evidence and model assumptions for each region are discussed in the following subsections. The 
goal is to find the magnetic field configuration within current observational bounds which results in a maximal 
photon- ALP mixing. In this way, a lower limit on can be derived. 

A. Magnetic field of the Milky Way 

The regular component of the B field of the Milky Way will be described with the analytical GMF model presented 
in [33]. The model consists of three components: a disk, a halo, and a so-called X component; and it predicts a 
field strength of the order (9(/iG). The model parameters were determined with a minimization utilizing the data 
of the WMAP7 polarized synchrotron emission maps and Faraday rotation measurements of extragalactic sources. 
Compared to previous models (e.g. Ref. [34 ), a relatively large field strength and extent is predicted for the halo 
and X component which leads to a comparatively large conversion probability in certain regions in the sky. The 
turbulent component of the GMF is neglected here, since the typical coherence length is of the order of 0(10 pc) 
which is far smaller than the oscillation length of photon- ALP-conversions. For each extragalactic VHE-7-ray source, 
the conversion probability is evaluated along the line of sight where it is assumed that the GMF is constant and 
homogeneous on a length scale of 100 pc. It was checked that smaller values for the domain length do not alter the 
results. Moreover, the density of the thermal electron plasma is calculated with the NE2001 code (in accordance with 
Ref. [33^) which predicts densities of the order of 10~^ cm~^ [35]. For further details on the conversion in the GMF, 
see Ref. [31 , where the same GMF model was utilized to compute the photon- ALP conversions. 

B. Intergalactic magnetic field 

In contrast to the GMF, little is known about the intergalactic magnetic field. From the observational side, only 
upper limits exist on the field strength, which constrain the IGMF at 2: = to a few 10~^ G [36]; and, e.g, the authors 
of Ref. [37 find ^iqmf = ^igmf(^ = 0) ^ 6 x 10~^ G for a coherence length of AjQj^p = 50Mpc using Faraday 
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rotation measurements of quasars. However, large scale structure formation with magnetic field amplification and 
cosmic ray deflection simulations suggest smaller values no larger than ^iqmf = 2 x 10~^^ G [38 or ^iqmf ^ 10~^^ G 
in voids [39 . The morphology of the IGMF is not known either and the most simple assumption is a domain-like 
structure which is also adopted here. The field strength is constant in each cell and only grows with cosmic expansion, 
i.e. Bigmf{z) = ^igmf(1 + but the orientation changes randomly from one cell to another. Furthermore, the 
domain length is given by AjQj^p. As shown in Ref. [27 , adopting a Kolmogorov-type turbulence spectrum instead 
of the simple domain structure has negligible effects on the results. In principle, the same procedure is followed here 
as presented in Ref. [26 , with the exception that the assumption of a the strong mixing is dropped. 

A scan over a logarithmic grid with 100 x 100 pixels in the {^igmf ^ ^igmf) space is performed to determine the 
most optimistic magnetic field setup. For each grid point, the photon survival probability is calculated with Eq. ( p!2| ) 
for 5000 realizations of the orientation of ^igmf for a fixed source distance z = 0.536, energy E = 0.574 TeV (this 
combination of z and E corresponds to an optical depth of r = 4 with the KD model), an ALP mass rria = 0.1 neV, 
and two different values of the coupling. Only the conversion in the intergalactic magnetic field with absorption due 
to the EBL of the KD model is taken into account. The impact of the photon- ALP conversions is quantified with the 
boost factor S, defined by 

B = P^^/exp(-r), (13) 

where P^^ is the median of the distribution of photon survival probabilities [40 . The result is shown for two different 
values of the photon- ALP-coupling in the top row of Fig. [1] As one would naively expect, for a large coupling of 
9aj = 5 X 10~^^GeV~^ (top-right panel) the boost factor increases with increasing -Sjqmf increasing AjQj^p up 
to a value of 10^'^ ^ 2.5. For even higher values, the boost factor starts to decrease again and shows an oscillatory 
behavior. This feature was already observed in Ref. [26^: if the conversion probability becomes too high, the photon 
fraction in the beam is large at all times; but at the same time, the photon flux is attenuated by the interaction with 
the EBL. As a consequence, B declines, and one is tempted to choose the values of AjQj^p and ^iqmf froi^i within the 
0.4 contour. The situation changes, however, if is decreased (top-left panel) by more than an order of magnitude 
to lO-^^GeV"^ The entire region of S > is shifted towards higher values in the {^igmf ^ -^igmf) plane. Without 
any a priori assumption about values of the ALP mass and coupling, it is thus advisable to select the maximum 
values of AjQj^p and ^iqmf ^^^^ allowed by observations and it is settled for AjQj^p= 50Mpc and ^iqmf — 5nG. 
For the thermal electron density in the intergalactic medium, a typical value of rigi, igm = 10~^cm~^ is adopted [41]. 



C. Intracluster magnetic fields 

In contrast to intergalactic magnetic fields, the existence of intracluster magnetic fields is well established. Syn- 
chrotron emission of the intracluster medium together with Faraday rotation measurements at radio frequencies have 
led to the common picture that turbulent magnetic fields of the order of (9(/iG) fill the cluster volume (see e.g. Refs. 
[42] [43] for reviews and typical values of the model parameters used below). The turbulence is usually described with 
a Kolmogorov-type spectrum, or with the simpler cell-like structure which is again used here (as in Ref. [31 ). There 
is evidence that the magnetic field strength follows the radial profile of the thermal electron distribution nei, icm in 
the cluster, 

BlCMF{r) = B^cMF (^ei,iCM(r)/n^i^icM)'^ . (14) 

with typical values 0.5 < r] < 1 and central magnetic fields up to ~ 10 /iG in the most massive clusters. The thermal 
electron density is described by 

^el,ICM(r) = n^ijCM (1 +^Acore)"^^^^ , (15) 

with characteristic values of /3 = 2/3 and rcore = 200 kpc. The coherence length is usually assumed to be comparable 
to galactic scales of the order of 10 kpc. As before, a grid scan over the (Aj^j^p, ^icmf) V^^^^ is performed in order to 
determine the parameters that maximize the photon- ALP conversions. A cluster with a radius of 2 Mpc is assumed 
together with r] = 0.5. Instead of the boost factor, the fraction of ALPs Pa^ in the final state [i.e., pfinai = Pss = 
(0, 0, 1)] is shown in the bottom row of Fig. [ijfor an initially unpolarized pure photon beam, Pa^ = Tt{ps3U p^npoi^^) • 
Again, 5000 ^icMF-field realizations are simulated, and the median Pa^ is computed. The more ALPs leave the 
cluster the stronger the effect will be on the VHE spectra because more ALPs can convert back into photons in the 
GMF and enhance the observed flux. The bottom-left panel of Fig. [l] clearly shows that more ALPs are produced 
for stronger magnetic fields and longer coherence lengths for a photon- ALP coupling strength of = 
Interestingly, the conversion probability is constant for constant values of Bq x -y/Lcoh (red dashed line), and Pa-y 
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Figure 1. Parameter space scan in the (Ac, Bo) plane. Top row: Photon- ALP-conversion in the IGMF. The color map displays 
the boost factor of the median of all simulated BioMF-field realizations; see Eq. (13). The adopted values for the coupling (Left 
column: large couplings; Right column: small couplings) are displayed in the figure together with the critical energy above 
which the conversion occurs in the strong mixing regime. Bottom row: Conversion in the ICMF. In the left panel, the color 
coding shows the fraction of the initial photon beam that is converted to ALPs (median over all realizations). The median of 
the conversion probability is constant for constant values of Bq x V^Z^oh, as indicated by the red dashed line. The bottom-right 
panel displays the dependence of Pa^ on the coupling ga-f for different values of Bq x V^Z^oh • In this panel, 68 % of all i^-field 
realizations for each Bq x VLcoh value fall into the corresponding shaded regions. See text for further details. 



increases quadratically with growing coupling strength until the maximum probability of ~ 1/3 is reached (bottom- 
right panel of Fig. [l]). Not surprisingly, the maximum is reached for smaller couplings for larger values of Bq x ^/L^. 
Thus, for an optimistic scenario, a central ICMF value of Bq = 10 /iG with a coherence length of 10 kpc is chosen 
and it is assumed that the VHE-emitting AGN is located at the center of a galaxy cluster. The core thermal electron 
density is taken to be n^^ icm ~ 10~^ cm~^. 



D. Magnetic field scenarios 



Now that the most optimistic values for the different magnetic fields are identified, four scenarios are presented for 
which the effect of photon- ALP-oscillations on VHE-7-ray spectra will be investigated. In all four configurations, the 
conversion in the GMF is included. 

1. In the first scenario, called general source hereafter, no specific environment is assumed for the ALP pro- 
duction and only the conversion in the GMF is included. Instead, an initial beam polarization pinit — 
1/3 diag(e~'^, e"'^, 1) is considered. This situation corresponds to a maximal mixing in some turbulent magnetic 
field inside or around the source and a subsequent attenuation of the photon fraction of the beam. In this 
general scheme, one is not forced to apply some sort of averaging over the many possible orientations of the 
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Table L Model parameters for the different magnetic field scenarios. In frameworks including the conversion inside galaxy 
clusters, the beam is assumed to travel the distance rduster through the volume filled with a B-field. In the optimistic 
ICM scenario, the B field varies as in Eq. (14). All AGN are assumed to be located at the center of a cluster. In the 



fiducial case, the magnetic field and thermal electron density are assumed to be constant throughout the cluster volume. Only 
AGN listed in Table 1 of Ref. ^SV are assumed to lie within galaxy clusters. See text for further details. 



Name 


IGMF 

^IGMF ^IGMF ^el,IGM 

(nG) (Mpc) (xlO-^cm-3) 


ICM 

^ICMF '^ICMF ^cluster ^el, ICM ^core Tj 

(/xG) (kpc) (Mpc) (xlO-^cm-2) (kpc) 


general source 


Only conversion in GMF, but pmit = l/3diag(e ^, e ^,1) 


optimistic IGMF 
optimistic ICM 
fiducial 


5 50 1 
0.01 10 1 


10 10 2 10 200 0.5 
1 10 2/3 1 



random magnetic field. 

2. In a second configuration, named optimistic ICM^ it is optimistically assumed that all VHE 7-ray-emitting AGN 
are located at the center of galaxy clusters of a 2 Mpc radius. The magnetic field changes over the distance 



from the cluster core as in Eq. (14) with a central magnetic field of ^icmf — 10 A^G and a coherence length of 
AjQj^p = 10 kpc. Any conversion in the intergalactic magnetic field is neglected, as well as any attenuation of the 
photon fiux by local radiation fields inside the galaxy cluster. Upon exit of the galaxy cluster, the photon beam 
will be attenuated by the interaction with the EBL whereas the ALP fraction propagates unhampered over the 
entire distance to the Milky Way. In the GMF, ALPs and photons can again convert into each other. Apart 
from the ICMF that drops radially with distance from the cluster core, this is the same setup as investigated in 
Ref. [31]. 

3. Thirdly, it will be assumed that no AGNs are affected by the photon- ALP-conversion inside a galaxy cluster 
but, on the other hand, the intergalactic field will be taken to its most optimistic values, i.e., ^iqmf — 5nG 
and AjQi^p = 50 Mpc. This setup is labeled optimistic IGMFd^nd is basically the same as that considered in, 
e.g., Ref. [26] apart from the complete energy-dependent treatment applied here. 

4. Finally, a set of more conservative model parameters is chosen to study both the conversion in the IGMF and 
ICMF. The parameters are conservative in the sense that they are not as close to the observational bounds as 
in the optimistic scenarios introduced above. Only the AGN listed in Tab. 1 of [31] are assumed to be located 
inside a galaxy cluster. As their position relative to the cluster core in unknown, a constant ICMF of 1 /iG is 
assumed. Furthermore, a value of rduster = 2/3 Mpc is adopted as the distance that photons propagate through 
the intra-cluster medium [44]. The value of the IGMF is motivated by simulations of large scale structure 
formations [38l [39] . This framework will be called fiducial. 

All scenarios are analyzed with two EBL models, namely the model of Ref. [9 (henceforth the FRV model) and the 
lower limit prediction of the KD model. The optical depth of the former is additionally scaled by a factor of ~ 1.3, as 
suggested by recent studies of VHE-7-ray spectra [45 . These two EBL models more or less bracket the range of the 
EBL density allowed by lower limits from galaxy number counts mentioned in Sec. |l| and upper limits derived from 
VHE-7-ray spectra [13]. Moreover, it was shown in Ref. [46] that these two models result in the highest significance 
of the pair production anomaly; and it can be expected that comparably small photon-ALP couplings are able to 
reduce this tension significantly and, thus, to derive conservative lower limits on the photon-ALP coupling. The 
different scenarios and their corresponding model parameters are summarized in Table [l| The photon- ALP-conversion 
inside the source is not explicitly taken into account here, but a possible contribution is accounted for in the general 
source scenario. 



IV. PROBING THE OPACITY WITH VHE GAMMA-RAY SPECTRA 



With the framework to calculate the photon survival probability P^^ introduced in the previous sections, the 
observed VHE-7-ray spectra are corrected for absorption in the presence of ALPs. As it is not assumed that the 
photon-ALP-conversions occur in the strong mixing regime, P^^ can show a strong oscillatory behavior. Therefore, 
an observed spectral point with a flux ^f^^ over an energy bin /S.Ei with central energy Ei is corrected with an average 
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transfer function, 

{P^^h = -^[ dEP^^{E), (16) 
SO that the absorption corrected flux <l>^ is obtained by 

= (p^^)-i $f ^ (17) 

In practice, the photon survival probability is calculated for 40 energies for each source and linearly interpolated in 
\ogiQ{E) and log2o(^77)- This has been cross-checked for one {rria^ga-f) pair with 100 energies and the results are 
found to be compatible if only 40 energies are used. 

The same technique as put forward in Appendix B of [18 is used here to quantify the significance of the pair 
production anomaly in the presence of ALPs. Each spectrum j with data points that correspond to r{z^Ei) ^ 2, i.e. 
the optical thick regime, is fitted with an analytical function fj{E). A list with all considered spectra that fulfill this 
criterion is shown in Table [ll} The function fj{E) is either a power law, or, in case the fit probability is pfn < 0.05, a 
logarithmic parabola, 

f(F\= iME/Eo)-^, Pflt > 0.05, 

' \iVo(E/So)-<^+''°'"<^/^°», otherwise, ^ ' 

with a flux and energy normahzation A^o and Eq, respectively; a power-law index F; and a curvature /3c. For each 
data point in the optical thick regime, the residual is calculated, 

which is normalized to the statistical measurement uncertainty on the fiux (68% confidence). Under the 
hypothesis that P^^ gives a correct prediction of the opacity of the Universe to VHE 7-rays, the residuals in the 
optical thick regime should follow a Gaussian distribution with zero mean. This conjecture is checked with the t test, 
for which the variable 

t= , (20) 

with X the mean and the variance of the residual distribution that contains data points, follows a t distribution 
from which the significance (one-sided confidence interval) can be calculated. This method to quantify the accordance 
between model and data has several advantages. Firstly, the functions to parametrize the spectra do not depend on 
any particular blazar emission model, as no constraints on the photon index T nor on the curvature Pc are made 
during the fit. Most spectra are adequately described by these functions, as shown in the Appendix. Secondly, no 
extrapolation from the optical thin to the optical thick regime is required, and the statistical uncertainties of the 
measurement enter the significance test self-consistently. 

Without the contribution of ALPs, one finds for the spectra listed in Table [ll] a significance of 5.1 x 10~^ ~ 4.4 <j 
for the KD model and 2.1 x 10~^ ^ 3.5 a for the FRV model scaled by 1.3 that the models do not describe the data. 
It might come as a surprise that the scaled FRV model gives a lower significance than the minimal attenuation model 
of Ref. [10 . The reason for this is that more data points migrate into the optical thick regime as the EBL density 
increases. This leads to an overall distribution closer to a zero mean and shows the limitation of the method: as long 
as the overall fit to all spectra is acceptable, the entire residual distribution must scatter around zero (see also Ref. 
146). 

In the following discussion, ALPs are included in the correction of the observed spectrum. For this purpose, the 
transition probability for all four scenarios and the two different EBL models is calculated separately for each source 
listed in Table [ll| This is necessary because each AGN has a different redshift (important for the attenuation) and a 
different position in the sky (influencing the conversion in the GMF). Furthermore, P^^ is computed over a grid of 
equally spaced values in the (log;Lo(^a)7 logio(^a7)) space. For the coupling constant, the range 10~^^ GeV~^ ^ ^ 
10-1° GeV"i is chosen for ah magnetic field frameworks. The upper bound is motivated by the bound set by the CAST 
experiment of < 8.8 x 10" GeV"^ [68] while for the lower bound the contribution of ALPs is expected to become 
negligible. On the other hand, the range of the tested ALP masses differs in the different scenarios. It is determined 
by the critical energy given in Eq. ([6| that should span an interval that includes the minimum and maximum energies 
of the VHE spectrum sample in Table [HI namely 0.08 TeV and 23.1 TeV. The different magnetic fields and thermal 
electron densities result in different mass ranges. For the general source- and optimistic ICM configurations, a mass 
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Table II. List of VHE-7-ray spectra included in the analysis. The table shows the redshift of the source, the lACT experiment 
that measured it, the energy range covered by the spectrum and the number of data points in the optical thick regime for the 
optical depth given by the KD model and by the scaled version of the FRV model. 



j 


Source 


Redshift 


Experiment 


Energy range 
(TeV) (r 


Nr^2 

= 1 X tkd) 


Nr^2 

(r = 1.3 X TFRv) 


Reference 


1 


Mrk421 


0.031 


HEGRA 


0.82 - 13.59 





1 


m 


2 


Mrk421 


0.031 


HEGRA 


0.82 - 13.59 





1 


m 


3 


Mrk421 


0.031 


H.E.S.S. 


1.12 - 17.44 





2 


m 


4 


Mrk421 


0.031 


H.E.S.S. 


1.75 - 23.10 


1 


4 


m 


r 



Mrk 50 


0.034 




U.OD - zl.4o 


1 


Q 


[50| 


6 


lES 1950+650 


0.048 


HEGRA 


1.59 - 10.00 





1 


m 


7 


lES 1950+650 


0.048 


HEGRA 


1.52 - 10.94 





1 


ED 


8 


PKS2155-30r 


[ 0.116 


H.E.S.S. 


0.23 - 2.28 





2 


EH 


9 


PKS 2155-304^ 


[ 0.116 


H.E.S.S. 


0.23 - 3.11 





3 


m 


10 


PKS 2155-304^ 


[ 0.116 


H.E.S.S. 


0.22 - 4.72 





6 




11 


PKS 2155-30^^ 


" 0.116 


H.E.S.S. 


0.25 - 3.20 





2 


EH 


12 


RGB J0710+591 0.125 


VERITAS 


0.42 - 3.65 





2 


m 


13 


H 1426+428 


0.13 


HEGRA, 
CAT, WHIPPLE 


0.25 - 10.12 


2 


5 


EH 


14 


lES 0229-200 


0.140 


H.E.S.S. 


0.60 - 11.45 


3 




M 


15 


H 2356-309 


0.165 


H.E.S.S. 


0.18 - 0.92 





1 


m 


16 


H 2356-309 


0.165 


H.E.S.S. 


0.22 - 0.91 





1 


EH 


17 


H 2356-309 


0.165 


H.E.S.S. 


0.23 - 1.71 





1 


m 


18 


lES 1218+304 


0.182 


VERITAS 


0.19 - 1.48 





3 


m 


19 


lES 1101-23^ 


0.186 


H.E.S.S. 


0.18 - 2.92 


3 


7 


m 


20 


lES 0347-121 


0.188 


H.E.S.S. 


0.30 - 3.03 


2 


4 


m 


21 


RBS041^ 


0.190 


VERITAS 


0.23 - 0.61 





1 


m 


22 


lES 0414+OOc 


I 0.287 


H.E.S.S. 


0.17- 1.13 


2 


3 


m 


23 


lES 0414+OOc 


0.287 


VERITAS 


0.23 - 0.61 





1 


m 


24 


PKS 1222+21 


0.432 


MAGIC 


0.08 - 0.35 





1 


ESI 


25 


3C 27^ 


0.536 


MAGIC 


0.15 - 0.35 


1 


1 


m 


26 


3C 279 


0.536 


MAGIC 


0.08 - 0.48 


1 


2 


115. 



^ Assumed to be located in a galaxy cluster in the fiducial scenario [see 1311 Tab. 1]. 
b Not included in [I8]. 



range of 1 neV ^ ^ 10^ neV is chosen; whereas for the optimistic IGMF setup, the smaller values of ^igmf and 
the ambient density lead to a shift in the mass to 10~^-^ neV ^ rria ^ 10^"^ neV. In the combined fiducial scenario, it 
is settled for the intermediate range 10~^-^ neV ^ ^ 10^-^ neV. A resolution of the grid of 32 x 32 = 1024 points 
is selected in the particular ranges of (log;Lo(^o), logio(^o7))- 

A complication is introduced by the random magnetic fields in the scenarios apart from the general source case. 
Since the exact orientation of the IGMF and ICMF in each domain is unknown, a large number A^sim of simulated 
random realizations is required. Here, A^sim will be set to 5000, and, therefore, for each {rria^gaj) pair one ends up with 
5000 values for the significance level of the t test, pt. One solution would be to compute the median (or mean) of the 
transfer function and afterwards calculate pt. However, in the averaging process all information on the pt distribution 
is lost, and it is unclear if this certain value is statistically suitable to deduce a lower limit on ^^7- Instead, the pt 
distribution is used to determine the pt value for which 95 % of all B-fieid realizations result in a worse compatibility 
of the particular framework with the data (i.e., those realizations that result in a smaller pt value). This particular 
significance is henceforth denoted as ^95. In summary, for each scenario, one now has one ^95 value for each grid 
point in the {rria^ga-f) space. The lower limit on ga^ is then defined as the contour line for which ^95 = 0.01. In this 
way, for {rria^gaj) values below this contour line, 95% of all 5-field realizations result in compatibility of less than 
1 % with the data. For the two EBL models used here, this corresponds to a decrease of the significance of the pair 
production anomaly by a factor of 5.1 x 10~^ (KD model) and 2.1 x 10~^ (scaled FRV model). 
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Figure 2. Significance map for the photon- ALP conversion in the {ma,ga-f) plane. Smaller values (brighter regions) indicate 
less accordance between the model and the data. Upper panel: pt values for the general source- scenario, shown as — log^Q(pt). 
Lower panel: ^95 values for the optimistic ICM case. For each pixel, 5 000 realizations of the random magnetic field are 
simulated and P95 is determined from the resulting 5000 pt values (cf. Sec. IV). In the left column, the attenuation due to the 
interaction of VHE 7-rays with the EBL is given by the KD model, whereas in the right column the FRV model is utilized. 
The maps are smoothed using a bilinear interpolation between the single (ma^Qa-i) pixels. 



V. RESULTS 



The results for the significance test introduced in the previous section are presented for each of the four scenarios 
developed in Sec. |III[ The upper panels of of Fig. |2] show the pt values for the general source configuration for 
the KD model (top-right panel) and the FRV model (top- left panel). In this scenario, no random magnetic field 
is involved, and thus there is only one pt value for each pixel. The color coding and the contour lines show the 
— \og^Q{pt) values and larger values of pt [smaller values of — log^LoCi^t)]' which represent a higher probability that 
the corresponding t value is the result of a statistical fluctuation; i.e., a higher probability that the transfer function 
is in accordance with the data. Clearly, the pt values increase with an increasing photon- ALP-coupling. The lower 
limit on {pt = 0.01 corresponds to the — \ogiQ{pt) = 2 contour line) is ~ 7.8 x 10~^^ GeV~^ for the KD model- 
and ~ 1.4 X 10~^^ GeV~^ for the FRV model, respectively, in the regime where the mixing becomes independent of 
the ALP mass at rria ^ 15neV. This mass marks the onset of the strong mixing regime (SMR) for all spectra in the 
environment of the Milky Way. For higher masses, the critical energy increases, and so does the number of spectral 
points outside the SMR. Higher couplings of are necessary to compensate this effect and to retain a low level of the 
significance of the pair production anomaly. Above ^ 250 neV, the tested coupling does not lead to a reduction 
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of the tension between the model and data in comparison to the no-ALPs case. These observations are vahd for both 
EBL models. 

A similar overall behavior is found in the optimistic ICM case (Fig. [2j bottom-right panel: KD model; bottom- left 
panel: FRV model). The color code now displays the ^95 values for the 5000 simulated realizations of the random 
ICMF in each pixel. Apart from the overall trend, peculiar regions are visible for the contour lines. In certain mass 
ranges, the lower limit contour for extends down to almost Qa^ = 10~^^ GeV~^ using the FRV model. These 
features are caused by the oscillatory behavior of P^^ outside the SMR which affects the low-energy data points in the 
spectra. These data points usually have the best count statistics, smallest error bars, and the strongest influence on 
the overall spectral fit. The oscillations in the transfer function can lead to a correction that is strong in one energy 
bin but small in the adjacent bin. As a result, the spectral fit is altered and leads to residuals in the optical thick 
regime that are closer to zero for certain {rua^gaj) pairs. 

Thus, it is expected that these features will change if more VHE spectra are included in a future analysis. Further- 
more, the oscillations of the transfer function lead to a poor fit quality for the spectra with the best overall statistics 
(Mrk421 [49], Mrk501 [50], and PKS 2155-304 [54 ) and to a smah overah fit probability (see the Appendix). This 
will lead to a broadening of the residual distribution and a possible overestimation of the ^95 values closer to 1. These 
features should not be taken as a preferred parameter region for ALPs to explain the opacity of the Universe. 

In the optimistic IGMF-scenaiio with the KD model, the only significant improvement over the no-ALP case is 
actually outside the SMR, as can be seen from Fig. [s] (top-right panel). Note that the mass range in which the 
transition to the SMR occurs has now shifted to lower masses due to the smaller IGMF and ambient electron density 
compared to the intracluster case. With the attenuation of the FRV model, the optimistic parameter choices for 
^IGMF ^IGMF ^^^^ ^ lower limit on ~ as low as ~ 3 x 10-^^ GeV"^ (top-left panel of Fig. m. 

The bottom row of Fig. |3] displays the results for the more conservative parameter choice of the fiducial-h dLmewoik. 
In this scenario, one cannot strictly speak about a lower limit on as neither the values of the magnetic fields nor the 
values for the coherence lengths are set to their observationally allowed upper limits. The {rria^ga-f) pairs that result 
in P95 ^0.01 can thus rather be seen as a preferred region in the parameter space if one tries to explain the opacity of 
the Universe with photon-ALP-conversions. One has to keep in mind, though, that the majority of simulated 5-field 
realization results in smaller pt values. Not surprisingly, one can conclude that the photon- ALP conversion in the 
IGMF is negligible, since the ^95 =0.01 contour line does not extend to lower values of ga^ at ^ 1 neV, as observed 
in the optimistic /GMF-case. Compared to the optimistic /CM-case, the lower limit contour line has shifted towards 
higher values in because a smaller number of AGNs is assumed to be located inside galaxy clusters. However, it 
has to be underlined that ALPs are still able to improve the accordance of the model with the data significantly. 

From Figs. [2] and [3] it is obvious that in the KD model, higher values of the photon- ALP-coupling are necessary to 
reduce the tension between model and data below the threshold of J995 =0.01 (more stringent lower limits) than in 
the scaled FRV model. The reason for this is twofold: On the one hand, without ALPs, the absorption correction in 
the FRV model is larger for high optical depths, which leads to higher residuals in some spectra. Lower photon- ALP- 
couplings suffice in these cases to reduce the residuals. On the other hand, the significance of the pair production 
anomaly is lower in the scaled FRV model to begin with (cf. Sec. jlVj ). Demanding the same decrease of significance 
as in the FRV model- without ALPs to the lower limit value (2.1 x 10~^ to 0.01) in the KD model-case results in 
a significance value of ~ 2.4 x 10~^, close to the pt = 10~^ contour line. Especially in the optimistic IGMF and 
fiducial scenarios this line is in good agreement with the pt = 10~^ contour line in the FRV model case. 

VIIE-7-ray spectra are subject to systematic uncertainties which can also affect the significance test used here. 
The authors of Ref. identified several sources of uncertainties in the quantification of the significance of the pair 
production anomaly such as a selection bias of VHE sources, the uncertainty of the overall energy scale of lACTs, 
and spillover effects in the highest energy bins due to the limited energy resolution of lACTs (the reader is referred 
to Ref. [H] for a detailed discussion). Including these effects leads in general to a reduction of the significance. The 
strongest reduction, to 2.6 <j, was found if the last energy bins of all spectra were excluded from the analysis and the 
energy points were simultaneously scaled by -15 % in energy (a conservative choice, as it was shown that a scaling 
of the order of 5 % is in better agreement with a cross correlation between lACTs and the Fermi-L AT [69 ). This 
certainly poses a lower limit on the significance, as it seems unlikely that all VHE spectra are influenced by these 
systematics in the same way. Nevertheless, a lower limit of ^95 = 0.01 with the inclusion of ALPs would not help 
to significantly improve the accordance between model and data in this case of a marginal indication. If ALPs were 
required at all, higher photon- ALP couplings would be necessary. However, the goal here is to set lower limits on gaj 
if the pair production anomaly is not explained by invoking all systematic uncertainties on the VHE observations at 
once. 

Figure [4] compares the lower limits derived here with current observational upper limits, regions of theoretical 
interest, and sensitivities of planned experiments. Only the lower limits for the FRV model are shown, since they all 
lie below the limits derived with the KD model. The lower limits clearly extend below the stringent upper limits from 
the CAST experiment [68^ (dark shaded region). In the optimistic IGMF case, they also lie below the upper limit 



12 




logio(m„ / neV) logio("^a / neV) 

Figure 3. Significance maps for the photon-ALP-conversion. Same as the bottom row in Fig. [2] but for the optimistic 
/CMF- scenario (top row) and fiducial-case (bottom row). 

derived from the nonobservation of prompt 7-rays from the supernova SN 1987a (gray shaded region) [70l[7l]. These 
7-rays would be the result of ALPs reconverted in the GMF that are produced in the supernova explosion [72] . The 
dotted-dashed lines show theoretical upper limits on calculated from magnetic white dwarfs (mWDs) [73] . Photon- 
ALP conversions lead to a linear polarization Pl of the photon beam [32 , and by treating the current observations of 
mWDs as a limit, i.e., Pl ^ 5%, one can derive a limit on the photon- ALP-coupling. The different lines correspond 
to different values of the magnetic field strength of the mWDs and different values for the limit on P^. Although the 
magnetic field and ambient density in the vicinity of mWDs are very different from the scenarios considered here, the 
mWD considerations turn out to be sensitive in the same (772^,^07) region as the VHE observations. Nevertheless, 
the limits use a P-field model inferred from one single mWD. If they are confirmed with future observations, they will 
strongly constrain the parameter space for ALPs that can potentially decrease the opacity of the Universe for VHE 
7-rays. 

The lower limits of the optimistic scenarios extend into the preferred region for the ALP parameters to explain 
the white dwarf (WD) cooling problem. It is difficult with current theoretical models to satisfactorily reproduce the 
observed WD luminosity function. Including the production of ALPs, on the other hand, with a mass and coupling 
within the light-gray-shaded band in Fig. [4] [74 serves as an additional cooling mechanism for WD and can reduce 
the tension between current model predictions and data [75 . This issue is, however, subject to ongoing discussion 

Interestingly, the ALP parameter space of the fiducial-scendino can be probed with planned experiments. The 
sensitivity forecasts for the improved Any Light Particle Search (ALPS II) |77l ^] and the International Axion 
Observatory (lAXO) [79 are displayed as a crosshatched and righthatched region, respectively, in Fig. [4j The lower 
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Figure 4. ALP parameter space with the lower limits on Qa-^ derived in this paper. The lower limits for the different scenarios 
are displayed as blue shaded regions, or in the case of the general so'izrce-scenario, as a dark blue solid line. They are only 
shown for the scaled FRV model, so that the optical depth is given by r = 1.3 x tfrv- For comparison, upper limits, hints for 
theoretical preferred regions, and sensitivities of future experiments are also plotted. See text for further details. 



limits derived here thus pose an additional physics case for these future experiments. 



VI. SUMMARY AND CONCLUSIONS 

In this article, VHE-7-ray observations have been used for the first time, to set lower limits on the photon- ALP cou- 
pling. Various magnetic field configurations have been considered, namely the magnetic field of the Milky Way, the 
intergalactic magnetic field, and the B field that pervades the intra-cluster medium. Under the assumption of domain- 
like random B fields, the photon- ALP-conversion probability for different field strengths and coherence lengths is 
investigated. As a result, hypothetical ALPs from a large region of the {rria^ga-f) parameter space would be able to 
reduce the opacity of the Universe for VHE-7-rays and to ease the tension between the predictions of common EBL 
models and data from lACTs that has been found in Ref. [18]. If field strengths and the coherence lengths of the 
magnetic fields are assumed to be close to their maximally observationally allowed values, lower limits on are 
obtained that reach down to ~ 10~^^GeV~^. Even lower values of ga^ could be obtained if, for instance, the 
GMF model were tuned to the most optimistic values allowed by the fitting errors calculated in Ref. [33 . Additional 
contributions to the GMF are also possible, such as a kiloparsec-scale magnetized wind [80 that could further enhance 
the conversion probability and reduce the lower limits. For more conservative values of the intervening B fields, the 
limits are of the order of ^ 2 x 10~^^GeV~^, close to the upper bounds of the CAST experiment and within the 
sensitivity estimates of future experiments such as ALPS II or lAXO. 



14 



Alternative mechanisms that are, in principle, able to increase the transparency of the Universe are also discussed 
in the literature. Active galactic nuclei could also be the source of cosmic-ray protons that produce secondary photons 
in the interaction with the EBL [81 . In such scenarios, the secondary photons are responsible for the high-energy 
end of VHE 7-ray spectra. 

The effect of the unknown EBL density on the lower limits of has been assessed with two different EBL 
models. At the time being, the sample of VHE spectra is dominated by sources with a redshift 0.1 < z < 0.2 and 
is thus most sensitive to changes in the EBL density at near-infrared wavelengths. As a consequence, certain EBL 
model realizations exist for which the pair production anomaly is less significant [46 , and higher values of Qa^ would 
be required to obtain a significant improvement over the situation without ALPs. One improvement would be to 
parametrize the EBL model independently (for instance, with splines, as done in, e.g., Ref. [13 ) and recalculate the 
significances in the presence of ALPs. This is left for future investigations. Firm conclusions can only be drawn with 
future direct observations of the EBL and VHE measurements in the optical thick regime of both distant sources 
at several hundreds of GeV and nearby sources at several tens of TeV, which will also enable further tests of ALP 
scenarios. Several such observations have already been announced [82 and will become more feasible with the next 
generation of air shower experiments such as the Cherenkov Telescope Array [83 , the High Altitude Water Cherenkov 
Experiment [84^, and the Hundred*i Square-km Cosmic ORigin Explorer [85j . 
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Appendix A: Fit qualities 

In this appendix, the quality of the fits of the analytical functions introduced in Eq. [18] is addressed. The fit 



statistics without the contribution of ALPs are listed in Table III For both EBL models, all spectra show a high fit 
probability. The only exceptions are the spectra of MrkSOl and one spectrum of 3C 279. The former spectrum is 
measured with high accuracy resulting in very small statistical errors. It is dominated by its systematical uncertainties 
(see gray band in Fig. 10 in [50 ) which are not included in the fit here. The spectrum of 3C 279 consist only of 
three data points, making a power law the only meaningful fitting function. The combined values translate into 
satisfactory overall fit probabilities of pfit = 0.160 and pfit = 0.303 for the KD model and FRV model, respectively. 

The fit qualities for the de-absorption of the observed VHE spectra with ALPs are poor in the transition to the 
strong mixing regime with p^t <C 1. The corresponding high values are dominated by a few spectra only (as in the 
no- ALPs case), namely those with high count statistics and consequently small error bars. The largest contributions 
come from the spectra of Mrk421 [49] and MrkSOl [50 , for the KD model, and, additionally, the spectrum of 
PKS 2155-304 [54] for the FRV model. Again, these spectra are dominated by their systematic uncertainties which 
are not included here. 

The reason for the large contribution of these spectra to the total values is the oscillatory behavior outside 
the SMR. As a result of the oscillations, the fit residuals also scatter strongly around zero and give rise to a low fit 
probability. If the spectra with high statistics are removed from the samples, the x^/d.o.f. values are close to one for 
ah {ma,ga-f) values. 

To summarize, the fits to most spectra in the {rria^ga-f) parameter space are acceptable. In the case of a small 
overall fit probability, the lower limits on are pushed towards lower values as the residual distribution broadens. 
Nevertheless, this bias mainly affects the region of ALP parameters outside the SMR and could be eliminated if the 
systematic uncertainties were included in the fit. 
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j 


oourcG 


rL/xpGrimGrit 


Fit 
functiorj^j 


(d.o.f.) 

T = 1 X 


Pfit 
tkd 


(d.o.f.) 

T = 1.3 X 


Pfit 
tfrv 


1 


Mrk421 


HEGRA 


LP 






4.10 (7) 


0.768 


2 


Mrk421 


HEGRA 


PL 






8.75 (8) 


0.364 


3 


Mrk421 


H.E.S.S. 


LP 






14.75 (10) 


0.142 


4 


Mrk421 


H.E.S.S. 


LP, Plj^ 


16.97 (11) 


0.109 


13.95 (12) 


0.304 


5 


Mrk 501 


HEGRA 


LP 


29.78 (14) 0.008 


38.55 (14) 


0.000 


6 


lES 1950+650 


HEGRA 


PL 






0.78 (3) 


0.854 


7 


lES 1950+650 


HEGRA 


PL 






7.89 (6) 


0.247 


8 


PKS 2155-304 


H.E.S.S. 


PL 






9.57 (7) 


0.215 


9 


PKS 2155-304 


H.E.S.S. 


PL 






8.34 (8) 


0.401 


10 


PKS 2155-304 


H.E.S.S. 


LP 






11.48 (14) 


0.648 


11 


PKS 2155-304 


H.E.S.S. 


PL 






5.46 (3) 


0.141 


12 


RGB J0710+591 


VERITAS 


PL 






2.00 (3) 


0.573 






HEGRA, 
CAT, WHIPPLE 


PT 


9.06 (10) 


0.526 




U.OUl 


14 


lES 0229-200 


H.E.S.S. 


PL 


3.17 (6) 


0.787 


3.86 (6) 


0.695 


15 


H 2356-309 


H.E.S.S. 


PL 






3.56 (6) 


0.735 


16 


H 2356-309 


H.E.S.S. 


PL 






4.20 (6) 


0.649 


17 


H 2356-309 


H.E.S.S. 


PL 






3.70 (6) 


0.717 


18 


lES 1218+304 


VERITAS 


PL 






2.33 (5) 


0.802 


19 


lES 1101-232 


H.E.S.S. 


PL 


5.70 (11) 


0.892 


6.83 (11) 


0.813 


20 


lES 0347-121 


H.E.S.S. 


PL 


2.59 (5) 


0.763 


2.30 (5) 


0.807 


21 


RBS0413 


VERITAS 


PL 






0.02 (2) 


0.988 


22 


lES 0414+009 


H.E.S.S. 


PL 


1.86 (4) 


0.761 


2.95 (4) 


0.567 


23 


lES 0414+009 


VERITAS 


PL 






0.51 (2) 


0.776 


24 


PKS 1222+21 


MAGIC 


PL 






0.24 (3) 


0.971 


25 


3C279 


MAGIC 


PL 


3.46 (1) 


0.063 


3.95 (1) 


0.047 


26 


3C279 


MAGIC 


PL 


3.64 (3) 


0.303 


4.48 (3) 


0.214 


Combined 


76.25 (65) 0.160 


173.88 (165) 0.303 



Table HI. List of fit qualities for all VHE-7-ray spectra if no ALPs are included in the de-absorption of the spectra. The table 
shows the values, the degrees of freedom (d.o.f.), and the resulting fit probabilities pfn for both EBL models used here. See 
Table [H] for the references of each spectrum. 

^ PL = power law, LP = logarithmic parabola. 

^ The spectrum is fitted with a logarithmic parabola in the KD model and with a power law in the FRV model, respectively. 
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